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Abstract 

We present various Lattice Boltzmann Models which reproduce the effects of 
rough walls, shear thinning and granular flow. We examine the boundary layers 
generated by the roughness of the walls. Shear thinning produces plug flow with a 
sharp density contrast at the boundaries. Density waves are spontaneously gener- 
ated when the viscosity has a nonlinear dependence on density which characterizes 
granular flow. 
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1. Introduction 

Many fluids in our daily life have rather complex rheological behavior. Pastes, 
suspensions, liquid crystals, dense polymers and granular media are usually non- 
Newtonian and can exhibit many flow anomalies, prominent ones being shear 
thinning or thickening and the spontaneous formation of density fluctuations in 
granular flow. Within the framework of classical fluid dynamics it is in general 
not simple to take into account these nonlinearities. Therefore it is of interest to 
look for alternative techniques to model the behavior of complex fluids. 

As one alternative to the direct solution of the equations of motion the so 
called Lattice Boltzmann Models (LBM) have been proposed^^'^J. These models 
are deflned on a lattice with velocity vectors that can only point into few discrete 
directions and all have the same length. This simpliflcation is somewhat com- 
pensated by the fact that on each site one has more real degrees of freedom (six 
on a triangular lattice) than in the classical numerical techniques allowing for the 
deflnition of a local shear or a local rotation. 

Two important questions concerning the LBM models can be asked: 1. How 
well do they reproduce solutions of the phenomenological equations of motion, like 
Navier Stokes? and 2. How well do they reproduce nature? The first question 
has been extensively addressed in the literature'"*^"^]. If certain assumptions are 
made on the length and time scales over which the variables can change the in- 
compressible Navier Stokes equation can be derived using the Chapman Enskog 
scheme. It is, however, known that straightforward simulations of the LBM can 
give inhomogeneous densities violating this incompressibility restriction. In this 
paper we will investigate these spatial and temporal density fiuctuations in more 
detail. In fact, we want to address mainly the second question: Can LB models 
describe real phenomena like shear thinning, density waves or the perturbations 
arising from the roughness of walls? 

For that purpose we will investigate the flow through a pipe along which the 
particles are accelerated through gravity. We want to see what happens if the 
walls of the pipe are rough and study constitutive laws that produce plug flow and 
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clogging. The typical experimental materials to which our investigations should 
apply, are suspensions with shear thinning in the case of plug flow, and granular 
media in the case of clogging. 

In the following section we describe the model and the various variants used 
in this paper. The next section is devoted to the effects of wall roughness. Section 
4 discusses models that give shear thinning and section 5 presents data for a model 
that spontaneously generates density waves. 

2. Description of the model 

We consider a triangular lattice, and on each site x we have six real variables 
Ni{x,t), i = 1,..6, representing (counted counter clockwise) the densities of the 
particles going in the direction i of the lattice. (For convenience we will in the 
following omit the site index x and denote by N'^ the value of the particle density 
after collision.) One updating of the system ( t — > t + 1 ) is given by two steps: 
(1.) The collision step at which the six Ni are updated at each site through 

Ni = N,+\{Ni-Nl'') (1) 

and (2.) the propagation step at which each Ni is shifted to the site of the nearest 
neighbor in direction i. Eq. (1) produces a relaxation towards the equilibrium 
densities N^'^ which is numerical stable provided the relaxation constant —2 < 
A < 0. The value of A sets the kinematic viscosity of the fluid. The equilibrium 
densities are given by 

Nf' = ^{l + 2u-Ci + A{u- c,f - 2^^) (2) 
where p is the mass density at site x 

i 

Cj the unity vector along direction i and u the velocity vector at site x defined 
through the momentum density per site 

pu = Y,c,N, . (4) 

i 

3 



The equilibrium distribution N^'^ given in Eq. (2), is chosen to give mass and 
momentum conservation in the coUision step. The flow wiU be forced into the 
direction of the gravity which is pointing paraUel to the walls of the pipe. For 
that purpose an additional step is added after the collision step which is defined by 
Nl' = Nl + ^Ci- (pg). Periodic boundary conditions are imposed in the direction of 
gravity in which the system has a length of Li. In the perpendicular direction one 
has walls separated by L2 lattice spacings. The lattice orientation is such that one 
of the lattice directions is parallel to the walls. At the beginning of the simulation 
the average density p is fixed. It is an important parameter of the model which 
because of mass conservation stays constant in time. We initialize the system by 
having the same values of the A^^ on each site and then let the system evolve to 
its steady state. In the case of the stable flows steady state is reached after 2000 
or 3000 time steps. In the case of the unstable flows that develop density waves, 
the simulations might take up to 20000 time steps to reach steady state. 

The sites lying on the walls of the system only have two directions a and 
b. Usually two different collision steps can be applied on these sitcs^^l, either the 
specular condition, i.e. — Nf^ and A^ = A^, or the bounce-back condition, i.e. 
A^ ^ = Na,b- In the propagation step for these sites the direction in which the Aj 
are shifted is inverted. We want to be able to implement walls that are not smooth 
but rugged, i.e. that have (quenched) disorder. For this purposes we introduce a 
mixed boundary condition defined through 

A^ = xNa + (1 - x)Nb and A^ = (1 - x)Na + xNb (5) 

where x = and y is a random variable chosen from a homogeneous distribution 
between and 1. Setting a = will give the pure bounce-back condition whereas 
CK = 00 corresponds to the pure specular reflection condition. 

The relaxation parameter A depends on the material properties including the 
kinematic and the bulk viscosities. Usually complex fluids are phenomenologically 
described by "constitutive laws" given e.g. by the functional dependence of the 
viscosities upon the shear velocity, the density or the pressure. We want to investi- 
gate the effect of rather typical nonlinear constitutive laws on the flow properties. 
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Since an exact relation between A and the material constants is not known we will 
lean on some approximative arguments^^'^'^l that predict a vanishing bulk viscos- 
ity. In that case one can relate A directly to the kinematic viscosity v through 
A = — 1(0.25 + 2v)~^. We will consider two cases: (1.) is a function of the local 
shear rate f and (2.) z/ is a function of the local density p. The detailed functional 
forms used here will be described in sections 4 and 5. 

Our calculations were performed on a Connection Machine CM-2 at GMD 
(Bonn) using 32-bit precision. The program needs less than one minute to make 
50 updates of a system of size 1024^. The program was also benchmarked on a 
CM-5 at I.P.G. in Parist^l 

3. The effects of rugged walls 

As already mentioned, it is well known that LB models produce inhomo- 
geneities in the density when used to simulate for instance flow through a pipe. 
In the middle of the pipe the density is higher than at the walls ^^^^ by a factor 
1/(1 —tt^) . This is seen in fig. 1 which shows the density in a cross section through 
the pipe. For o; = 0, i.e. the case of smooth walls the density profile has precisely 
the predicted shape of pwaii/i^ — u"^) as can be seen from the line showing the 
pressure ^^^^ p = {p/2){l — v?). In the case of rough walls, for which we have 
chosen in ck = 1 in fig. 1, the density has a minimum close to the walls. Also, the 
pressure has a minimum at the wall and a lower, but still constant value in the 
center. As expected there are some random fluctuations close to the walls. 

In fig. 2 we see the density variation along the center of the pipe. Since the 
values taken at different times coincide very well one is in the steady state. Clearly 
the randomness of the refiection properties of the wall still have some effect but 
the relative variation is of the order of 0.0001, i.e. extremely weak. 

The roughness at the walls therefore seems to be screened very efficiently. 
This is seen more clearly in fig. 3 which shows the entire density profile in the 
pipe. The boundary layer has a thickness of a few mean free paths (the mean free 
path in this context is the characteristic length over which a perturbation in the 
NiS will be damped and has typical length of 1/A) where the distribution of the 
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Ni^s is clearly different from that in the bulk of the material. In this sense it may 
be characterized as a Knudsen layer. 

4. Shear thinning and plug flow 

Shear thinning can be phenomenologically explained by a non-Newtonian con- 
stitutive law given by a decrease of the viscosity as a function of the shear rate. 
Within the context of the LBM the shear rate f can be defined through 

H^) = ^\^CiU\\{x + Ci)\ (6) 

i 

where u\\ is the projection of u into the direction of the pipe. We consider a 
constitutive law of the form u = ui for t < tq and u = U2 for r > tq and ui > U2. 

In fig. 4 we show the velocity profile in a cross section through the pipe. In the 
simulations the flow was initialized with a relatively strong forcing, = 5 x 10~^. 
During this initial phase the shearthinned regions at the walls appear, and the 
flow velocity increases to approximately its steady state value. Then the forcing 
was reduced by a factor lOtogr^SxlO"^ and the system allowed to reach steady 
state. 

We see that the profile is rather flat in a broad central region which ends at 
a sharp kink after which one finds a rather steep velocity gradient towards the 
walls where the fluid is in the thinned phase. This kind of behavior is usually 
called plug flow. It was checked that the flow is really in a steady state by per- 
forming longer runs. In a recent preprint^^^ a simulation of a similar LBM has 
been presented which also finds plug fiow by taking a constitutive law in which 
the viscosity decreases with the shear rate like a power law. This seems to indicate 
that the appearance of plug flow is rather independent on the detailed form of the 
constitutive law as long as z/ is a decreasing function of f. 

5. Density waves 

A salient feature of granular media is the spontaneous formation of density 
waves, similar in fact to traffic jams on highways. One possibility to explain the 
effect that generates these waves is to assume that the viscosity depends on density. 
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Within the kinetic gas theory of granular mediat^'^'^^1 the relation z/ oc (p — Pc)^^^ 
has been derived. Since the above relation imposes a maximum density pc it is 
rather difficult to implement it directly within the context of the LBM where the 
particles do not have an exclusive volume. We therefore chose a piecewise linear 
relation of the form v = Vmin if P ^ Pt and v = i'o + j{p — 'p) for p > pt (see fig. 5). 
p is the average density and the threshold density pt is chosen to make u a positive 
continuous function of the density. Fig. 6 shows results from simulations where 
Pt = 2.962 and the slope 7 = 6.25 corresponding to a minimum cut-off viscosity 

l^min = 0.01. 

In order to generate density waves we found it necessary to introduce a small 
perturbation producing a 0.3% relative density difference. This perturbation was 
performed by introducing a small amount of momentum on one line across the 
pipe, keeping the mass unchanged. In fig. 6 we see that this initially very weak 
perturbation dramatically builds up and develops into a density wave of over 10% 
density contrast. For a pipe of same width but half the length, i.e. a different 
aspect ratio the wave has a less pronounced profile. This dependence on the 
aspect ratio is not to be confounded with finite size effects. Our mean free path 
is typically one lattice spacing so that the strong finite size effects encountered 
in some lattice gas models'-*^^! should not be relevant here. The maximum flow 
velocity Umax at the later times is Umax = 0.039 and Umax = 0.048 for the channel 
lengths 256 and 512 respectively and same width 64. The forcing g = 3.33 x 10~^ 
is the same for both system lengths. For the present parameter values the initial 
perturbation relaxes, leaving a time independent density field, if the value of 7 
is less than 3.75. This effect can be understood qualitatively by observing that 
there are two competing mechanisms in the system: On one hand, the viscous 
relaxation of density perturbations will tend to smoothen density contrasts. On 
the other hand, the rather steep increase of viscosity with density combined with 
the presence of the walls will tend to increase the contrasts. A small increase 
in the density at the wall will give a local increase in viscosity and slow down 
the flow. Due to the inertia of the surrounding flow this in turn will lead to a 
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further increase in the density and so on. If this instabihty dominates the relaxing 
mechanism the density wave wiU form. 

By triggering the density wave by two spatially separated perturbations, 
rather than just a single one, we checked that the complex shape of the waves 
does not reflect the detailed way in which they were initiated. We also observe 
that there seems to be no characteristic wavelength: Fig. 6 shows that the waves 
have roughly the same shape on the scale of the channel length although the 
amplitude depends on the system size. 

Fig. 7 shows this amplitude as a function of time during 60,000 time steps. The 
insert shows the initial unstable phase leading to the rather drastic increase of the 
amplitude at the time 10,000. The first small increase in the density is due to the 
accelleration of the flow and can be understood from the velocity dependence in the 
pressure. The small jump at the time 2500 results from the perturbation. Before 
the instability is triggered at the time 10,000, small oscillations in the amplitude 
are observed. It was checked that the amplitude indeed has its' steady state value 
at the time 60,000 by running the simulations ten times longer. The complicated 
relaxation towards the fully developed density wave indicates that strong non- 
linear effects come into play rendering a linear stability analysis meaningless. It 
would be interesting to understand this behaviour further. 

In flg. 8 one can see the density wave propagating. The fronts are actually 
sharpest at the boundary and the left gradient which is less sharp than the right 
one has some weak spatial oscillations. The fact that the waves are of the order 
of the length of the pipe again shows that there is no characteristic length scale. 
The periodic boundary conditions seem crucial to reinforce the steady state. One 
therefore has the typical behaviour of a kinetic wave as also found in traflic jam 
modelst^^l 

6. Conclusion 

We have presented various versions of Lattice Boltzmann Models which can 
reproduce rather complex flow behavior. On one hand we investigated how the 
laminar flow screens the asperities arising from rugged walls by forming Knudsen 
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layers close to the wall. When a shear thinning constitutive law is introduced 
we find plug flow. Finally, when the viscosity is an increasing function of density 
we observe a range of parameters for which the material spontaneously produces 
density waves traveling upstream. These density waves are triggered by some 
perturbation that apparently is unstable, but the final shape of the wave is inde- 
pendent of the initial disturbance. 

Plug flow and density waves are common phenomena in non-Newtonian fluid 
dynamics and have been investigated recently in detail for granular media'^^]. 
It seems therefore that LB models can be a powerful tool to handle complex 
fluids numerically. This approach is, however, yet quite preliminary. One has 
to determine the physical parameters for which the proposed models do match a 
real experiment and then compare measured and simulated results. Work in this 
direction is in progress. 
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Figure Captions 

Fig. 1 Density as a function of the position Y across the channel for a = (squares) 
and a = 1 (x). Also 2p = p{l — tt^) is shown for a = (o) and a = 1 (+). 
Li = 256, La = 64, p = 3, ^ = 3.33 x 10"^ and u = 0.25. The figure shows 
the steady state profile after 7500 iteration steps. The maximum flow velocity 

"^max ~ 0.52. 

Fig. 2 Density in the center of the pipe as a function of the position X along the 
channel for cu = 1 and otherwise the same parameters as in fig. 1. The two 
lines correspond to two measurements 25 iteration steps apart. 

Fig. 3 Density contrast in the pipe for the same parameters as in fig. 1 and fig. 2. 
White denotes the lowest density and black the largest one. 

Fig. 4 Velocity as a function of the position, Y across the channel measured at steady 
state after 5000 iteration steps . The insert shows the viscosity's dependence 
on the local shearrate. In this simulation vi = 1.0, 1^2 = 0.1, tq = 10~^, 
p = 3.0, Li = 256 and L2 = 64. The forcing is ^ = 5 x 10"^. 

Fig. 5 The density dependence of the viscosity chosen in the simulations. 

Fig. 6 The density in the center of the channel as a function of the position X along 
the channel for pt = 2.962, p = 3.0, g = 3.33 x 10"^ and L2 = 64. The 
curve of crosses is for Li = 256 and 60,000 iteration steps after the initial 
perturbation. The other curves correspond to Li = 512 and 5000 (thick line), 
60,000 (full line) and 60,025 (dashed line) iterations after the perturbation 
was applied. The slope 7 = 6.25 and the minimum viscosity Uj^in = 0.01. 

Fig. 7 The amplitude, i.e. difference between largest and smallest density, along the 
center of the pipe as a function of time measured in units of 100 iteration 
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steps for Li = 256 and otherwise the same parameters as in fig. 5. The insert 
is a blow-up of the behavior at early times. 
Fig. 8 Density contrast in the pipe for the same parameters as in fig. 5 and fig. 6 
after 60,000 (upper) and 60,025 (lower) iteration steps. White denotes the 
lowest and black the largest density. 
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